function HPZ_MME_Estimation (data,obs_num,treatment,flag,beta_flag,zeros_flag,aggregation)
% The function calculates the utility functions of the GARP satisfying 
% subjects in CFGK (2007) for a given bundle using NLLS.
% data is a matrix with six columns (We assume that the endowment is fixed at 1):
% The first column is the subject ID.
% The second column is the observation number - 50 observations per subject
% The third column is the quantity of good 1 chosen by the subject.
% The fourth column is the quantity of good 2 chosen by the subject.
% The fifth column is the price of good 1. 
% The sixth column is the price of good 2. 
% subjects is a vector that lists the subjects IDs.
% subject_num is the number of subjects.
% obs_num is the number of observations per subject.
% treatmentis the number of treatment in CFGK (2007).
% flag is the type of utility function chosen for the prizes:
% 1 - CRRA.
% 2 - CARA.
% The function returns a matrix with subject_num rows and nine columns: 
% The first column is the estimated beta using the max measure.
% The second column is the estimated additional parameter (rho if CRRA, 
% A if CARA) using the max measure.
% The third column is the value of the function at the optimum using the max measure.
% The fourth column is the estimated beta using the average measure.
% The fifth column is the estimated additional parameter (rho if CRRA, 
% A if CARA) using the average measure.
% The sixth column is the value of the function at the optimum using the average measure.
% The seventh column is the estimated beta using the sum of squares measure.
% The eighth column is the estimated additional parameter (rho if CRRA, 
% A if CARA) using the sum of squares measure.
% The ninth column is the value of the function at the optimum using the sum of squares measure.

%% This file keep track of all results produced by the fminsearchbnd
global fid_MME_final_report
%% "True" value of this boolean flag restricts beta to be equal to zero
global beta_zero_flag
%% "True" value of this boolean flag restricts rho to be equal to zero
global rho_zero_flag
%% The max number of convergance points MME
global MME_min_counter;
% Max number of initial points for the fminsearchbnd optimization procedure
global max_starting_points
% now we only use max 100 iterations as our upper bound
max_starting_points = 100;
%% TO BE DONE
% Time (minutes) that user want to spend on running the code on his machine
global MME_max_time_estimation

if MME_max_time_estimation == Inf
    MME_time = MME_max_time_estimation;
else
    % the entire time that user wants to invest on estimation
    MME_time = MME_max_time_estimation * 60; 
end

% keep track of the entire estimation process
t_total = tic;

% choose proper algorithm for the fminsearchbnd optimization procedure
options = optimset('Algorithm','interior-point','Display','off');

% subject_data is a matrix with obs_num rows and 4 columns.
% each row is one choice of the subject.
% The first column is the quantity of good 1 chosen by the subject.
% The second column is the quantity of good 2 chosen by the subject.
% The third column is the price of good 1. 
% The fourth column is the price of good 2. 
observations = data(1:obs_num,3:6);

if zeros_flag==2 % Choi et al. (2007) correction
    
    % Choi et al. (2007) correction is applied for corner choices
    subject_data = HPZ_No_Corners (observations,obs_num,1);
    
else
    
    subject_data = observations;

end

%% Handle the corner solutions when CRRA is chosen as the utility function 
restricted_rho = false;
% a flag that restrict rho to be in the range of [0,1] for subjects with
% cornor solutions where CRRA is chosen as the utility function
if ((any(observations(:,1) == 0) > 0 || any(observations(:,2) == 0) > 0) && flag == 1)
    % check whether subject has any corner solution(s) (meaning x1 or x2 = 0)
    % among their 50 observations, where the chosen utility function is CRRA
    restricted_rho = true; 
end

% rho = 0 case 
zero_rho_initial = false;
if (flag == 1 || flag == 2)
    zero_rho_initial = true; 
end

%% Setting initial point for the fminsearchbnd optimization procedure
% The first two initial points are the point chosen by CFGK (2007) and the
% rest are choosen points randomly.
% initial_points = HPZ_Initial_Points (1,restricted_rho);
initial_points = HPZ_Initial_Points (1,restricted_rho,zero_rho_initial);

% The endowments are normalized to 1
endowment=ones(obs_num,1);

% Initializing the structures that hold the estimation results
% (this matrix keep track of all estimations of parameters)
results=zeros(max_starting_points,2);
   
% (keep track of all fvals)
criterion=zeros(max_starting_points,1);

%% ASK USER TO PROVIDE THIS NUMBER
% % number of optimal fvals
% MME_min_counter = 100;

% Beta Rho Fval : 3xmin_fvals
optimal_parameter_matrix = zeros(MME_min_counter, 3);

% initialization of threshold on the aggregated criterion value
fval_threshold = 1e-5;

% number of consecutive equal fvals
equal_fval_counter = 0;

% define the waitbar
h_wb = waitbar(0,strcat('Recovering Subject-',num2str(data(1,1)), ' Preferences. Please wait...'), 'name', 'Estimation is running...');
new_bar_val = 0;

% This array keep track of time that is spent on each covergence point
time_optim_conv_pnt = zeros(1, MME_min_counter);

% keep track of accumulative time for all initial points
time_init_accum = 0;
time_init_pt = 0;

for j=1:max_starting_points
    
    if MME_time == Inf
        %% Estimation withought time limit
    
        if (equal_fval_counter < MME_min_counter)
            

            %% Updating the waitbar
            if equal_fval_counter > 0 
                % update the value of the progress bar
                new_bar_val = equal_fval_counter / MME_min_counter;
            end
            % computations take place here
            waitbar((j / max_starting_points) + new_bar_val);

            % calculate estimation time for each convergance point
            t_optim = tic;
            
            if beta_zero_flag == true && restricted_rho == false
                % if beta restricted to be equal to 0, and there is no
                % constraint on the rho's value, then:

                % set beta to be 0, as we restricted it to be so.
                results(j,1) = 0;
                % the optimization process only search for the optimal value
                % of the rho, given beta is equal to 0.
                % Note: in this case values that can be taken by rho are all
                % positive values within the range of (0,inf).
                [results(j,2),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper_Extreme_Param(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,2), 0, [], options);
                
            elseif beta_zero_flag == true && restricted_rho == true
                % if beta restricted to be equal to 0, and there is a
                % constraint on the rho's value, then:

                % set beta to be 0, as we restricted it to be so.
                results(j,1) = 0;
                % the optimization process only search for the optimal value
                % of the rho, given beta is equal to 0.
                % Note: in this case values that can be taken by rho are all
                % positive values within the range of (0,1).
                [results(j,2),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper_Extreme_Param(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,2), 0, 1, options);

            elseif rho_zero_flag == true && beta_flag==1
                % if rho restricted to be equal to 0, and beta can take negative
                % values, then:

                % set rho to be 0, as we restricted it to be so.
                results(j,2) = 0;
                % the optimization process only search for the optimal value
                % of the beta, given rho is equal to 0.
                % Note: in this case values that can be taken by beta are all
                % values within the range of [0,inf).
                [results(j,1),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper_Extreme_Param(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,1), 0, [], options);

            elseif rho_zero_flag == true && beta_flag==2
                % if rho restricted to be equal to 0, and and there is no
                % constraint on the beta's value, then:

                % set rho to be 0, as we restricted it to be so.
                results(j,2) = 0;
                % the optimization process only search for the optimal value
                % of the beta, given rho is equal to 0.
                % Note: in this case values that can be taken by beta are all
                % positive values within the range of [0,inf).
                [results(j,1),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper_Extreme_Param(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,1), 0, [], options);

            elseif beta_flag==2 && restricted_rho == false
                [results(j,:),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,:),[0 0], [], options);
            
            elseif beta_flag==2 && restricted_rho == true
                % if the subject has a corner solution(s) where CRRA is chosen 
                % as the utility function, then restrict rho<1 and beta can be
                % any number positive number [0, inf]
                [results(j,:),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,:),[0 0], [inf 1], options);

            elseif beta_flag==1 && restricted_rho == false
                [results(j,:),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,:),[-1 0], [], options);
                
            elseif beta_flag==1 && restricted_rho == true
                % if the subject has a corner solution(s) where CRRA is chosen 
                % as the utility function, then restrict rho<1 and beta can be
                % a negative number [-1, inf]
                [results(j,:),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,:),[-1 0], [inf 1], options);
            end

            % initialize the temporary minimum of fvals by using the objective 
            % value of the first initial point
            if (j == 1)
                fval_temp_min = criterion(j);
                % initialize equal_fval_counter to 1
                equal_fval_counter = equal_fval_counter + 1;
                % update the estimation running time for the convergence point
                time_optim_conv_pnt(equal_fval_counter) = toc(t_optim);

                optimal_parameter_matrix(equal_fval_counter, 1:2) = results(j,:);
                optimal_parameter_matrix(equal_fval_counter,   3) =  criterion(j);
            end

            % check if the resulting criterion is equal to computed minimum
            % criterion or not AND number of minimal criterion are not greater 
            % than min_fvals (=3) AND j > 1
            if ( abs(fval_temp_min - criterion(j)) <= fval_threshold && equal_fval_counter < MME_min_counter && j > 1)
                % update the convergence counter
                equal_fval_counter = equal_fval_counter + 1;

                % update the estimation running time for the convergence point
                time_optim_conv_pnt(equal_fval_counter) = toc(t_optim);

                % store new minimum in the matrix
                optimal_parameter_matrix(equal_fval_counter, 1:2) = results(j,:);
                optimal_parameter_matrix(equal_fval_counter,   3) =  criterion(j);

                % update the minimum when the threshold is not achieved
                if(criterion(j) < fval_temp_min)
                    % update the temporary minimum of fvals
                    fval_temp_min = criterion(j);

                    % check the previous results to see if they are within the
                    % new acceptance range for criterion

                    % number of violations
                    mismatch_counter = 0;
                    % copy the current number of convergance points
                    temp_equal_fval_counter = equal_fval_counter;

                    % among all current convergance points:
                    for j_ind = 1 : equal_fval_counter
                        % check if they all are within new acceptance range
                        % if not:
                        if (optimal_parameter_matrix(j_ind,   3) > (fval_temp_min + fval_threshold))
                            % replace that element with inf
                            optimal_parameter_matrix(j_ind,   :) = Inf;
                            % update the value of temporary counter for
                            % convergance points
                            temp_equal_fval_counter = temp_equal_fval_counter - 1;
                            % update the number of violations
                            mismatch_counter = mismatch_counter + 1;
                        end
                        % Do nothing
                    end

                    % re-arrange the results vector:
                    % start index
                    start_ind = 1;
                    % end index
                    end_ind = equal_fval_counter;
                    % temporary value for swaping values
                    temp_val = zeros(1,3);

                    % Go through the result vector as long as there is still
                    % violation
                    while mismatch_counter > 0 && end_ind > 0 && start_ind < equal_fval_counter
                        % push Infs to the end of the array
                        if (optimal_parameter_matrix(start_ind,   3) == Inf)
                            % swap the Inf element with the last element
                            temp_val = optimal_parameter_matrix(start_ind,   :);
                            optimal_parameter_matrix(start_ind,   :) = optimal_parameter_matrix(end_ind,   :);
                            optimal_parameter_matrix(end_ind,   :) = temp_val;
                            % update mismatch counter value
                            mismatch_counter = mismatch_counter - 1;
                            % update pointers
                            end_ind  = end_ind - 1;

                        end
                        % move pointer forward
                        start_ind = start_ind + 1;
                    end

                    % update the equal_fval_counter
                    equal_fval_counter = temp_equal_fval_counter;

                end


            % check if we find a new minimum AND number of minimal criterion
            % are not greater than min_fvals (=3) AND j > 1 AND criterion(j) <
            % fval_temp_min 
            elseif (criterion(j) < fval_temp_min && ...
                    abs(fval_temp_min - criterion(j)) > fval_threshold && ...
                    equal_fval_counter < MME_min_counter && ...
                    j > 1)
                % update minimum
                fval_temp_min = criterion(j);
                % update counter
                equal_fval_counter = 0;
                % update result matrix
                optimal_parameter_matrix = zeros(MME_min_counter, 3);

                % store new minimum in the matrix
                equal_fval_counter = equal_fval_counter + 1;
                optimal_parameter_matrix(equal_fval_counter, 1:2) = results(j,:);
                optimal_parameter_matrix(equal_fval_counter,   3) =  criterion(j);

            end

            % if optimization is terminated
            if equal_fval_counter == MME_min_counter || j == max_starting_points
                % store the total number of used initial points
                total_init_points = j;
            end

        end
    
    else
        %% Time limited estimation
           
        %% Updating the waitbar
        % computations take place here
        waitbar(time_init_pt/MME_time);
        
        % calculate estimation time for each convergance point
        t_optim = tic;
        
        % calculate estimation time for each initial point
        t_init = tic;
        
        if j <= 3 || ((j>=3 && j<=10) && (MME_time > (time_init_accum + time_init_pt)))

            if beta_zero_flag == true
                % if beta restricted to be equal to 0, and there is no
                % constraint on the rho's value, then:

                % set beta to be 0, as we restricted it to be so.
                results(j,1) = 0;
                % the optimization process only search for the optimal value
                % of the rho, given beta is equal to 0.
                % Note: in this case values that can be taken by rho are all
                % positive values within the range of (0,inf).
                [results(j,2),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper_Extreme_Param(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,2), 0, [], options);

            elseif rho_zero_flag == true && beta_flag==1
                % if rho restricted to be equal to 0, and beta can take negative
                % values, then:

                % set rho to be 0, as we restricted it to be so.
                results(j,2) = 0;
                % the optimization process only search for the optimal value
                % of the beta, given rho is equal to 0.
                % Note: in this case values that can be taken by beta are all
                % values within the range of [0,inf).
                [results(j,1),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper_Extreme_Param(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,1), 0, [], options);

            elseif rho_zero_flag == true && beta_flag==2
                % if rho restricted to be equal to 0, and and there is no
                % constraint on the beta's value, then:

                % set rho to be 0, as we restricted it to be so.
                results(j,2) = 0;
                % the optimization process only search for the optimal value
                % of the beta, given rho is equal to 0.
                % Note: in this case values that can be taken by beta are all
                % positive values within the range of [0,inf).
                [results(j,1),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper_Extreme_Param(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,1), 0, [], options);

            elseif beta_flag==2
                [results(j,:),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,:),[0 0], [], options);

            elseif beta_flag==1
                [results(j,:),criterion(j)]=fminsearchbnd(@(param) HPZ_MME_Helper(param,endowment,subject_data,treatment,flag,aggregation),initial_points(j,:),[-1 0], [], options);
            end

            % update the estimation running time for the convergence point
            time_optim_conv_pnt(j) = toc(t_optim);

            optimal_parameter_matrix(j, 1:2) = results(j,:);
            optimal_parameter_matrix(j,   3) =  criterion(j);

            % measure the time for each initial point
            time_init_pt = toc(t_init);
            % measure the accumulative time time for all initial points so far
            time_init_accum = time_init_accum + time_init_pt;
            % store the total number of used initial points
            total_init_points = j;
        end
    
    end
        
end
    
% update the waitbar
remaindar_multplier = 0.20;
new_bar_val = 1 - (remaindar_multplier * 1);
waitbar(new_bar_val);

%% TO Be Done
%% choose distinct results
% % using threshold on optimal parameter sets
% 
% % This variable keeps track of the number of distinc optimal parameter sets
% distinct_result_counter = 0;
% 
% % These two thresholds define which optimal parameter sets are distinct
% diff_dinstinc_threshold = 5e-1; % 0.5
% ratio_dinstinc_threshold = 5e-2; % 0.05
% 
% % For all optimal results compute the ratio and diffrence in order to find
% % distinct results
% for i = 1 : MME_min_counter
%     repmat(optimal_parameter_matrix(:, 1), size(optimal_parameter_matrix, 1), 1)
% end


%% compute all results
% [~,best_index]=min(criterion);
% output(1,1)=results(best_index,1);
% output(1,2)=results(best_index,2);
% disp(['The estimated parameters are beta = ',num2str(output(1,1)),' and roh\A = ',num2str(output(1,2)),]);
% disp('Conluding by evaluating other criteria for these estimated parameters');
% output(1,3)=HPZ_Second_Criterion(output(1,1:2),observations,treatment,flag,zeros_flag,2,1);
% output(1,4)=HPZ_Second_Criterion(output(1,1:2),observations,treatment,flag,zeros_flag,1,1);
% [output(1,5),output(1,6),output(1,7)]=HPZ_MME_Criterion(output(1,1:2),endowment,observations,treatment,flag);
%% REMOVED
% output(1,8)=HPZ_Budget_Lines_Test (subject_data,treatment,flag,output(1,1:2));

if MME_time == Inf
    
    final_output = zeros(MME_min_counter, 7);
    for i = 1 : MME_min_counter
        % [~,best_index]=min(criterion);
        final_output(i,1)=optimal_parameter_matrix(i, 1);
        final_output(i,2)=optimal_parameter_matrix(i, 2);

        % NLLS Criterion using Choi et al. (2007) metric
        final_output(i,3)=HPZ_Second_Criterion(final_output(i,1:2),observations,treatment,flag,zeros_flag,2,1);
        % NLLS Criterion using Euclidean metric
        final_output(i,4)=HPZ_Second_Criterion(final_output(i,1:2),observations,treatment,flag,zeros_flag,1,1);
        % MME Criterion (Max Waste, Mean Waste, Sum of Squares Wastes)
        [final_output(i,5),final_output(i,6),final_output(i,7)]=HPZ_MME_Criterion(final_output(i,1:2),endowment,observations,treatment,flag);

        %% finalize the output file for distinct results
        fprintf(fid_MME_final_report, '%f,%f,%f,%f,%f,%f,%f,%f,%s\n', ...
            final_output(i,1), final_output(i,2), ...
            final_output(i,3), final_output(i,4), ...
            final_output(i,5), final_output(i,6), final_output(i,7), ...
            time_optim_conv_pnt(i), num2str(data(1,1)));

        % update the waitbar
        waitbar(new_bar_val + ((i/MME_min_counter)*remaindar_multplier));
    end
    
else
    %% time limited estimation
    final_output = zeros(total_init_points, 7);
    for i = 1 : total_init_points
        % [~,best_index]=min(criterion);
        final_output(i,1)=optimal_parameter_matrix(i, 1);
        final_output(i,2)=optimal_parameter_matrix(i, 2);

        % NLLS Criterion using Choi et al. (2007) metric
        final_output(i,3)=HPZ_Second_Criterion(final_output(i,1:2),observations,treatment,flag,zeros_flag,2,1);
        % NLLS Criterion using Euclidean metric
        final_output(i,4)=HPZ_Second_Criterion(final_output(i,1:2),observations,treatment,flag,zeros_flag,1,1);
        % MME Criterion (Max Waste, Mean Waste, Sum of Squares Wastes)
        [final_output(i,5),final_output(i,6),final_output(i,7)]=HPZ_MME_Criterion(final_output(i,1:2),endowment,observations,treatment,flag);

        %% finalize the output file for distinct results
        fprintf(fid_MME_final_report, '%f,%f,%f,%f,%f,%f,%f,%f,%s\n', ...
            final_output(i,1), final_output(i,2), ...
            final_output(i,3), final_output(i,4), ...
            final_output(i,5), final_output(i,6), final_output(i,7), ...
            time_optim_conv_pnt(i), num2str(data(1,1)));

        % update the waitbar
        waitbar(new_bar_val + ((i/total_init_points)*remaindar_multplier));
    end
end

% keep track of the entire estimation process
time_total = toc(t_total);

% % % % %% Add an empty line to the file
% % % % fprintf(fid_MME_final_report, '%s,%s,%s,%s,%s,%s,%s,%s,%s\n', ...
% % % %                               '','', ...
% % % %                               '','', ...
% % % %                               '','','', ...
% % % %                               '','');
% % % % if aggregation == 1
% % % %     aggregation_str = 'Max Waste';
% % % % elseif aggregation == 2
% % % %     aggregation_str = 'Mean Waste';
% % % % elseif aggregation == 3
% % % %     aggregation_str = 'Sum of Squares Wastes';
% % % % end
% % % % 
% % % % %% Add Total time and total number of used initial points to the report
% % % % fprintf(fid_MME_final_report, '%s,%s,%s,%s,%s,%s,%s,%s,%s\n', ...
% % % %                               strcat('Subject-', num2str(data(1,1))), ...
% % % %                               'Total_Estimation_Time(sec)',num2str(time_total), ...
% % % %                               'Total_InitialPoints_Used',num2mstr(total_init_points), ...
% % % %                               '',strcat('Money Metric Method - ', aggregation_str), ...
% % % %                               '','');
% % % % 
% % % % %% Add an empty line to the file
% % % % fprintf(fid_MME_final_report, '%s,%s,%s,%s,%s,%s,%s,%s,%s\n', ...
% % % %                               '','', ...
% % % %                               '','', ...
% % % %                               '','','', ...
% % % %                               '','');
% % % % %% Add an empty line to the file
% % % % fprintf(fid_MME_final_report, '%s,%s,%s,%s,%s,%s,%s,%s,%s\n', ...
% % % %                               '','', ...
% % % %                               '','', ...
% % % %                               '','','', ...
% % % %                               '','');

% close the waitbar
close(h_wb);

end